#include "yocto/z++/gzstream.hpp"
#include "yocto/exception.hpp"
#include "yocto/sequence/vector.hpp"
#include "yocto/ios/ocstream.hpp"
#include "yocto/math/types.hpp"

#include <iostream>

using namespace yocto;
using namespace math;

typedef double Real;


static inline bool load_frame( double &stamp, vector<Real> &v, ios::istream &fp, double &tc )
{
    char C = 0;
    if( fp.query(C) )
    {
        fp.store(C);
        stamp          = fp.read<double>();
        const size_t n = fp.read<uint32_t>();
        if( n <= 0 ) throw exception("invalid #data");
        v.make(n,0.0);
        fp.load( &v[1], n * sizeof(Real) );
        tc = fp.read<double>();
        return true;
    }
    else {
        return false;
    }
}

void compare_db()
{
    
    
    vector<Real>   norm;
    vector<Real>   chem;
    ios::igzstream norm_in( "norm-db.gz" );
    ios::igzstream chem_in( "chem-db.gz" );
    ios::ocstream  fp(  "delta.dat",   false );
    ios::ocstream  fp2( "speedup.dat", false );
                       
    double norm_stamp, chem_stamp;
    
    //--------------------------------------------------------------------------
    //-- load normal frame
    //--------------------------------------------------------------------------
    double tc_norm,tc_chem;
    double sum = 0;
    size_t num = 0;
    while( load_frame( norm_stamp, norm, norm_in, tc_norm ) )
    {
        //----------------------------------------------------------------------
        //-- load corresponding chem frame
        //----------------------------------------------------------------------
        if( !load_frame( chem_stamp, chem, chem_in, tc_chem ) )
        {
            throw exception("Mismatch #frames");
        }
        
        if( chem.size() != norm.size() )
            throw exception("Mismatch sizes");
        
        //-- compute all
        Real dmax = Fabs( chem[1] - norm[1] )/( Fabs(chem[1]) + Fabs(norm[1]) );
        for( size_t i=1; i <= chem.size(); ++i )
        {
            const Real tmp = Fabs( chem[i] - norm[i] )/(Fabs(chem[i]) + Fabs(norm[i]));
            if( tmp > dmax )
                dmax = tmp;
        }
        
        fp("%g %g\n", norm_stamp, dmax );
       
        if( tc_chem > 0 )
        {
            const double ratio = tc_norm/tc_chem;
            sum += ratio;
            ++num;
            const double ave = sum/num;
            fp2("%g %g %g\n", norm_stamp, ratio, ave );
        }
        
    }
        
}
